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Abstract 



The computation of free energy differences through an exponential weighting of out-of- 
equiUbrium paths (known as the Jarzynski equaUty [15, 16j) is often used for transitions 
between states described by an external parameter A in the Hamiltonian. We present 
here an extension to transitions between states defined by different values of some reac- 
tion coordinate, using a projected Brownian dynamics. In contrast with other approaches 
(see e.g. [22]), we use a projection rather than a constraining potential to let the con- 
straints associated with the reaction coordinate evolve. We show how to use the Lagrange 
multipHers associated with these constraints to compute the work associated with a given 
trajectory. Appropriate discretizations are proposed. Some numerical results demonstrate 
the applicability of the method for the computation of free energy difference profiles. 

Keywords: free energy, mean force, constrained dynamics, sampling techniques, Jarzyn- 
ski equality, Feynman-Kac formula. 

The free energy of a system is a quantity of paramount importance in statistical physics. 
It is of the form 



where (5 = 1/(/cbT') (T denotes the temperature and /cb the Boltzmann constant) and Z is 
the partition function 



of the Boltzmann (or Gibbs) measure exp(— /3y)d/i. In this expression, the function V = V{q) 
is the potential energy of the system (denoting by q the position vector) and ^ is a reference 
positive measure with support S. The space S is the configuration space of the system. We 
will consider here that S is a submanifold of M^^, but all the results extend to the case 
when S is a submanifold of (the 3 A^- dimensional torus, which arises when using periodic 
boundary conditions). The statistics of the system are completely defined by {V,^). 

In most cases, iV, fi) is labeled using a d-dimensional parameter z (with d <C 3A^) which 
characterizes the system at some coarser level. The parameter z can be independent of the 
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current configuration of the system. In this case, only the expression of the potential V depends 
on the parameter, so that the associated switching has sometimes been called 'alchemical 
transition'. Some examples of such parameters are the intensity of an external magnetic field 
for a spin system, or the temperature for a simulated annealing process. However, it is often 
the case that the parameter z labels submanifolds of the configuration space, through level 
sets = {CiQ) = of some function ^. The function ^ is called a 'reaction coordinate'. In 
this case, /i (especially the support of /i) depends on z and is defined using the orthogonal 
projection from or T^'^ to (this will be made precise in Section II. ip . Standard 
examples of reaction coordinates are bond lengths or dihedral angles in a molecule. 

The absolute free energy ([1]) can be computed only for certain systems, such as ideal 
gases, or solids at low temperature (resorting to the phonon spectrum) [23]. However, in 
many applications, the quantity of interest is the free energy difference between an initial 
and a final state (characterized by two different values of the parameter z). The free energy 
difference profiles indeed give information about the relative stabilities of several species, as 
well as their transition kinetics. The free energy differences are much more amenable to 
compute than the absolute free energy. Classical techniques to this end fall within three main 
classes. The first one, dating back to Kirkwood [17] . is thermodynamic integration^ which 
mimics the quasi-static evolution of a system as a succession of equilibrium samplings, which 
amounts to an infinitely slow switching between the initial and final states. The second one, 
the free energy perturbation method, was introduced by Zwanzig [35]. It recasts free energy 
differences as a phase-space integral, so that usual sampling techniques can be employed. 
Notice also that there exist many refinements for those two classes of techniques, such as 
umbrella sampling [3l]. The last and most recent class of methods uses dynamics arising from 
a switching at a finite rate. This can be done using nonequilihrium dynamics (the so-called fast 
growth methods) with a suitable exponential reweighting, as introduced by Jarzynski in |15| . 
Notice that the thermodynamic integration and free energy perturbation methods can be seen 
respectively as the limits of infinitely slow and fast switching of nonequilibrium dynamics, at 
least formally. Instead of being imposed a priori, this switching may also arise as the result 
of an equilibrium sampling, using for example the Adaptive Biasing Force technique [TJ [12] or 
metadynamics [H]. In those cases, the system is progressively forced to leave regions where 
the sampling of the reaction coordinate has been completed. 

It is still a matter of debate which method is the most efficient. While some results 
show that fast growth methods can be competitive in some situations [H], other studies 
disagree [19]. The results of [19] indeed indicate that even with the use of efficient path 
sampling techniques (see also [SniEOlEl]), fast growth methods do not outperform conventional 
methods such as umbrella sampling or thermodynamic integration (at least in a number of 
typical cases). However, general conclusions about the efficiency of fast growth methods are 
still to be drawn, depending on the cases under consideration. We believe that there is room for 
improvements of this relatively new method {e.g. by optimizing the switching schedule |24|). 
Let us also mention that this method is straightforward to parallelize and naturally provides 
with a posteriori error bounds via the central limit theorem, since it involves many independent 
trajectories. 

Most methods to compute free energy differences are well suited to the alchemical transition 
setting, but do not straightforwardly extend to the reaction coordinate setting. This latter 
case is the focus throughout this article. In this case, the methods described above require 
to consider dynamics restricted to the submanifold S^. For computations using Hamiltonian 
dynamics, we refer for example to [Hill]. In the stochastic case, thermodynamic integration 
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in the reaction coordinate case using projected stochastic dynamics has recently been put on 
a firm grounding [HIE]- On the other hand, stochastic nonequihbrium dynamics d la Jarzynski 
in the reaction coordinate case was, to our knowledge, not studied mathematically. It is the 
aim of this paper to perform such a study and to present a methodology to compute free 
energy differences in this framework. 

Nonequihbrium computations of free energy differences in the reaction coordinate setting 
using stochastic dynamics have until now used soft constraints to switch between the initial 
state centered on the submanifold {(,{q) = zq} and the final state centered on {£,{q) = zi}. 
Steered molecular dynamics techniques use for example a penalty term K{^{q) — z)"^ in the 
energy of the system [22] (with K large) to 'softly' constraint the system to remain close to the 
submanifold {£,{q) — z = 0}, and varying the value z from to 1 in a finite time T . It is shown 
in [13] how to use such a biasing potential to exactly compute free energy differences (even 
for a finite K), which is of particular interest for experimental studies. Prom a computational 
viewpoint however, it is expected that large values of K require small integration time steps. 
Moreover, it is observed in practice that the statistical fiuctuations increase with larger K 
(see [22]) ■ Instead, we propose to replace the stiff constraining potential K{^{q) — z)"^ by a 
projection onto the submanifold {(,{q) — z = 0}. This situation is reminiscent of the case of 
molecular constraints, that can be enforced using a stiff penalty term, or more elegantly and 
often more efficiently, using some projection of the dynamics involving Lagrange multipliers. 
This is the spirit of the well known SHAKE algorithm |26| . 

We propose a nonequihbrium stochastic dynamics and an equality that allow to compute 
free energy differences between states defined by different values of a reaction coordinate. The 
dynamics relies on a projection onto the current submanifold at each time step, and we use 
the Lagrange multipliers associated with this projection to estimate the free energy difference. 
More precisely, we use the difference between these Lagrange multipliers and the external 
forcing term required for the finite time switching (see for example the discretization (I3ip ). 
The main results of the paper are the Feynman-Kac equality of Theorem 12.21 (which extends 
the proof of [13] to hard constraints), as well as the associated discretizations ([33]) and (f34l) . 

The method we propose forces the system to pass free energy barriers, and thus enables 
free energy difference computations for metastable systems. Of course the reliability of the 
algorithm crucially depends on the choice of the reaction coordinate, which represents the 
essential degrees of freedom. The reaction coordinate should be rich enough in order to 
adequately describe the configuration paths of the system from the initial state to the final 
state. The determination of the essential degrees of freedom of a system is a very important 
problem, which is not the focus of this work. Thus, in the following, we suppose that a "good" 
reaction coordinate is given, and we are interested in the computation of free energy differences 
associated with this reaction coordinate. 

Let us also notice that some recent refinements of nonequihbrium dynamics to compute free 
energy differences, especially path sampling techniques [34] and Interacting Particle Systems 
approaches [25] (which equilibrate the nonequihbrium dynamics through some birth/death 
process based on the current work), can be extended to the reaction coordinate setting using 
the techniques we present here. Moreover, we restrict ourselves to the so-called overdamped 
Langevin dynamics, but it is possible to extend these results to the usual Langevin dynamics 
(this is a work in progress). 

The paper is organized as follows. In Section [U the thermodynamic integration setting is 
outlined in the reaction coordinate case. Section [2] then extends the method to nonequihbrium 
dynamics. Adapted numerical schemes are proposed in Section [3], and some numerical results 
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assessing the correctness of the method are presented in Section [H For clarity, we present the 
method in the case of a one-dimensional reaction coordinate and postpone until Appendix El 
the proofs and the expressions for the multi-dimensional case. 



1 Equilibrium computation of free energy differences 

The aim of this section is to introduce the definitions of the free energy and the mean force, 
and to recall how thermodynamic integration is used to compute free energy differences. The 
computation of the mean force is based on projected stochastic differential equations (SDE). 
These SDEs will also be used for the discretization of Jarzynski equality in Section [2l This 
section mainly reviews results of [6]. 



1.1 Free energy and mean force 



In the following, we denote by C M the configuration space of the system when no 
parameter z is involved. The state of the system is characterized by the value of a reaction 
coordinate : — > [0, 1]. The function is supposed to be smooth and such that S/^{q) / 
for all q £ ^A. For a given value z G [0, 1], we denote by the submanifold 

= { g G A^, ^q) = ^ } (3) 

and we assume that Uzep i] <^ For each point g G S^, we also introduce the orthogonal 
projection operator P{q) onto the tangent space to at point q defined by: 

Ffa)=Id-^^(,). (4) 

where denotes the tensor product. The orthogonal projection operator on the normal space 
to T,z at point q is defined by P'^{q) = Id — P{q). 
The free energy is then defined as 

Fiz) = -p-hn{Z,), (5) 

with 

Zz = 1^ eM-PV)daE., (6) 



where for any submanifold S of M , cjs denotes the Lebesgue measure induced on S as a 
submanifold of M^^. The associated Boltzmann probability measure is 

d^l^^=Z-'exp{-(]V)daJ:^. (7) 

Remark 1.1 (On the definition of the free energy). Two comments are in order about for- 
mula 1^. First, this formula is valid up to an additive constant, which is not important when 
considering free energy differences. Second, the potential V in (S^ may be a potential different 
from the actual potential seen by the particles. More precisely, if the particles evolve in a 
potential V, the standard definition of the free energy in the physics and chemistry literature 
is 1^ with 

Zz = j exp(-/3F) 5g(g)-2, 
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where is a measure supported by and defined by: for all test functions 4>, 

This amounts to considering with V replaced by an effective potential V + /3~^ Iti\\/^\ 

(see Remark \A.l\ for the case of a multi- dimensional constraint). Since the results we present 
in this paper hold irrespective of the physical signification of the potential V , we may assume 
without loss of mathematical generality that the free energy is indeed given by IS^-^B^. Let 
us emphasize that, in practice, the cumbersome computation of the gradient of the additional 
term /3~^ln|V^| in the modified potential (which intervenes in the projected SDEs we use, 
see /[i)7|)-/ f^) or ^29 ) - ^^) ) can be avoided resorting to some finite differences, as explained 
in ^. 

Using the co-area formula (see (H2l) and Proposition lA.2l for a proof in the multi-dimensional 
case), it is possible to derive the following expression of the derivative of the free energy F 
with respect to z (the so-called mean force) (see [211 EI]): 



where 



(8) 



is the mean curvature vector field of the surface S^. The free energy can thus be expressed as 
an average with respect to //s^: 

F'(z) = ^ fiq)dfij:M, (10) 
where / is the local mean force defined by: 

/ = ^-(V^ + /3-'^)- (11) 

In next section, we will explain how it is possible to compute this average with respect to fi-£^ , 
without explicitly computing /, by using projected SDEs. This avoids in particular the com- 
putation of the mean curvature vector H which involves second-order derivatives of ^. 
The principle of thermodynamic integration is to recast the free energy difference 

AF{z) = F{z) - F{0) (12) 

between two reaction coordinates and z as an integral over the mean force: 

AF{z) = r F'{y) dy. (13) 







Therefore, in practice, thermodynamic integration computation of free- energy is as follows. 
First, the free energy difference AF(z) is estimated using quadrature formulae for the integral 
in (fT3]) . such as for example a Gauss-Lobatto scheme: 



K 



AF{z)o^J2'^iF'{y, 



i=0 
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where the points {yo, yi, . . . , yx} are in [0, z] and {loq, loi, . . . , ujk} are their associated weights. 
Second, the derivatives F'{yi) are computed as canonical averages over the submanifolds S^^. , 
using projected SDEs (see next section). 

To obtain a free-energy profile (and not only a free-energy difference for a fixed final state), 
it is possible to approximate the function AF{z) on the interval [0, 1] by a polynomial. This 
can be done for example by interpolating the derivative F' by splines, and integrating the 
resulting function (consistently with the normalization AF(0) = 0). 



1.2 Projected stochastic differential equations 

In this section, we explain how to compute the mean force F'{z) defined by ([8]) using projected 
SDEs, for a fixed parameter z. We consider the solution Qt to the following SDE: 

^ , (14) 

where Bt is the standard 3 A^- dimensional Brownian motion and o denotes the Stratonovich 
product. It is possible (see [6]) to check that /is^ is an invariant probability measure associated 
with the SDE (fT^ . Under suitable assumptions, which we assume in the rest of the section, on 
the potential V and the surface S^, the process Qt is ergodic with respect to //s^- Moreover, 
the SDE (fT^ can be rewritten in the following way: 



dQt = -VV{Qt) dt + VW^dBt + VaQt)dAt, (15) 

where Aj is a real valued process, which can be interpreted as the Lagrange multiplier asso- 
ciated with the constraint (,{Qt) = z (see the discretization in Section [3TT]) . This process can 
be decomposed into two parts: 

dAt = dKf + dK\. (16) 
The so-called martingal^ part A™ (whose fiuctuation is of order \/At over a timestep At) is 

dkT = -VW^^iQt) ■ dBt, (17) 

where • implicitly denotes the Ito product. The so-called bounded variation part A^ (whose 
fluctuation is of order At over a timestep At) is 

d^t = p^(Qt) • VV^(Qt) dt + /?-i^(Qi) . H{Qt) dt = f{Qt) dt, (18) 



/ being the local mean force deflned above by (fTTl) . Thus, since Qt is ergodic with respect to 
the mean force can be obtained as a mean over the Lagrange multiplier At: 

Proposition 1.2. The mean force is given by: 

F'{z) = hm 1 r dAt = hm 1 dA^ (19) 

T^oo 1 Jq T^oo 1 Jq 



^For our purposes, it is enough to think of a martingale as an Ito integral with respect to the Brownian 
motion {Bt)t>o- 
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Notice that the martingale part dA™, which has the largest fluctuations, has zero mean. 
In order to reduce the variance, it is thus numerically convenient to perform the mean over the 
bounded variation part dAl rather than over the whole Lagrange multiplier dAt (see Section [S]). 

We refer to [6] for a proof of Proposition 11.21 as well as for formulae involving higher 
dimensional reaction coordinates. Such ideas have been used for a long time in the framework 
of Hamiltonian dynamics (see [2n [27] ). 

The interest of Equation (fTO]) is that the SDE (fT5]l can be very naturally discretized 
as explained in Section 13.11 below. Then, the average over a discretized trajectory of the 
process At converges to F'{z). This is particularly convenient for numerical purposes since 
it does not ask for explicitly computing the local force /. For further details, we refer to [6] 
and to Section 13.11 In next section, we use these ideas for the computation of the free energy 
difference given through the Jarzynski equality. 

2 Nonequilibrium stochastic methods in the reaction coordi- 
nate case 

As opposed to quasistatic methods where the free energy difference between an initial state 
and a flnal state is expressed by (fT3]) . in nonequilibrium methods, the free energy difference 
is expressed using a Feynman-Kac average over nonequilibrium paths [13 dSj [25] 

AF(1) = F(l) - F(0) = hiE (^e-'^^(^)) , (20) 

where W(T) denotes the total work exerted along a nonequilibrium path {Qt, •z(i))te[o,T] > with 
2;(0) = and z{T) = 1. 

We wish here to extend the Feynman-Kac formula derived in [13] for a parameter z which 
appears only in the potential V , to the reaction coordinate case, where z labels submanifolds 
(deflned by Equation ([3])) of the state space. To this end, we need to make precise the evolution 
of the constraints. 

We consider a path z : [0, T] — > [0, 1] of values of the reaction coordinate ^, with 
z(0) = 0, and z{T) = 1. Recall that the associated family of submanifolds of admissible 
configurations is denoted by 

and that the associated Boltzmann probability measures are 

dMs,(t) = ^7(J)exp(-^y)das,„). 

We construct a diffusion {Qt)t(^[Q,T] so that Qt G ^^(t) for all t G [0,T] and {Qt)t(^[Q,T\ satisfies 
the following properties (see Section [2?n for a more rigorous formulation): 

• For all t G [0, T], Qt+dt is the orthogonal projection on T,z{t+dt) of the position obtained 
by the unconstrained displacement: Qt — W{Qt)dt + ^j2j3~^dBt. 

To each realization of this process, a work W{t) can be associated as 

W(t) = f f{Qs)z'{s)ds, 
Jo 
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where / is the local mean force defined above by ([TTI) . Then, we prove that the Feynman- 
Kac formula (QUI) holds for the free energy F associated with the reaction coordinate and 
defined by ([5|). Notice that, at least formally, in the limit of an infinitely slow switching 
from z{0) = to z(T) = 1, Formula (f20l) corresponds to the thermodynamic integration 
formula (fTSl) . Formula (f2Ql) enables the computation of free energy difference at arbitrary 
rates, through a correction consisting in a reweighting of the nonequilibrium paths. 

The rest of this section is organized as follows. In Section 12.11 we make precise the pro- 
cess Qt we consider. Then, in Section 12.21 we state the Feynman-Kac formula (f20]l for a 
one-dimensional reaction coordinate. We recall that the formulae for the general case in- 
volving higher dimensional reaction coordinates, as well as the main proofs, are presented in 
Appendix [Al 



2.1 The nonequilibrium projected stochastic dynamics 

The considered diffusion reads, in the Stratonovich setting: 

dQt = -P{Qt)yV{Qt)dt + ^/W'P{Qt)odBt + Vi{Qt)dM'^' 



At) ' ' ^^'^ 



With a view to the discretization of Qt, let us notice that Qt can be characterized by the 
following property: 

Proposition 2.1. The process Qt solution to (|21|) is the only ltd process satisfying for some 
real-valued adapted ltd process {J^t)te[o,T]- 



(0)' 



dQt = -VViQt)dt + ^/W'dBt + VC{Qt)dAt, 

aQt) = z{t). 

Moreover, the process {J^t)t£[{),T] &e decomposed as 

At = AJ^ + A[ + Af*, (22) 

with the martingale part 



dhT = -^j2p=^^{Qt)-dBt, 
the local force part (see i fTT]) for the definition of f) 

= • (^^(Q*) dt + r^H{Qt)) dt = f{Qt) dt, (23) 



and the external forcing (or switching) term 

* |ve(Q 



dAr' = dt. 
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The proof of Proposition 12.11 is easy and consists in computing d^{Qt) by Ito's calculus 
and identifying the bounded variation and the martingale parts of the stochastic processes. 

The difference with the projected stochastic differential equation (fT4l) considered in the 
thermodynamic integration setting is that the out-of-equilibrium evolution of the constraints z{t) 
creates a drift V^{Qt) dA^^^ along the reaction coordinate. This drift can be interpreted as 
an external forcing required for the switching to take place at a finite rate, and must be sub- 
tracted from the Lagrange multiplier At in order to obtain a correct expression for the work 
Wit) involved in the Feynman-Kac fluctuation equality (see Equations ([3TI1 and ((331) below). 
This correction is quantitatively important when the switching is not slow. 



2.2 The Feynman-Kac fluctuation equality 

Let us define the nonequilibrium work exerted on the diffusion (f2T]) by: 



mi)= t f{Qs)z'{s)ds, (24) 
Jo 

where / is the local mean force defined above by ([TT]) . In practice, the nonequilibrium 
work W{t) can be computed by using the local force part dAl (see (f23]l ). as in the ther- 
modynamic integration method (see (fT9]l ). Thus, the formula we use to compute yV{t) is 
rather: ^ 

Wit) = [ /(s)(iAi, (25) 
Jo 

since can be obtained by a natural numerical scheme (see Section E]), avoiding the cum- 
bersome computations of the mean curvature vector H in the expression of / (as already 
explained in Section [TTT]) . 

We can now state the generalization of the Jarzynski nonequilibrium equality to the case 
when the switching is parameterized by a reaction coordinate. 

Theorem 2.2 (Feynman-Kac fluctuation equality). For any test function ip and\/t G [0,T], 
it holds 

^ I ^d^E.,,=EMQ,)e-/^^W). 
In particular, we have the work fluctuation identity: Vt G [0,?"], 

AF{z{t)) = F{z{t)) - F{z{0)) = In (e fe'^^W)) . (26) 



As in the alchemical case [13| . the proof follows from a Feynman-Kac formula. The proof 
of this theorem is presented in the general multi-dimensional case in Appendix [A] (see Theo- 
rem lAlSl). 



3 Discretization of the dynamics 

The main interest of the above formulae (fT3]) - (fT9]) and (f25]) - (f26]l is that they admit natural 
time discretizations. The principle is to use a predictor-corrector scheme for the associated 
dynamics (fT^ and (PTI) . and to use the Lagrange multiplier Aj to compute the local mean 
force /. 
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Section 13.11 is mainly a review of the results of [6] and presents this idea in the context of 
thermodynamic integration. Then, we extend the method to the case of evolving constraints 
in Section [321 



3.1 Discretization of the projected difFusion 

For the projected SDE (fTSl) onto a submanifold = {S,{q) — z = 0}, two discretizations of the 
dynamics, extending the usual Euler-Maruyama scheme, are proposed in [S]. These numerical 
schemes for constrained Brownian dynamics are in the spirit of the so-called RATTLE p] and 
SHAKE [26] algorithms classical used for constrained Hamiltonian dynamics, and also related 
with the algorithms proposed in [32l [Tl [20]. 
The first one is: 

Qn+l = 

where AA^+i is such that ^(Qn+i) = z, 

where At is the time step and ?7" is a 3A^-dimensional standard Gaussian random vector. 
Notice that ([271) admits a natural variational interpretation, since Qn+i can be seen as the 
closest point on the submanifold to the predicted position Qn — ^ViQn) At + \/2Atp~^ Un- 
The real AA„+i is then the Lagrange multiplier associated with the constraint ^(Q„+i) = z. 
Another possible discretization of ([T5l) is 



Qn+i — Qn - yViQn) At + yj2At(5-^Un + AA„+i VC(Q„), 
where AA^+i is such that ^(Qn+i) = z. 

Although this scheme is not naturally associated with a variational principle, it may be more 
practical since its formulation is more explicit. Notice also that we use the same notation AA^ 
for the Lagrange multipliers for both ([27ll and ([28ll (and later for ([29ll and ([30])), since all the 
formulas we state in terms of AA„ are verified whatever the constrained dynamics. 

To solve Equation ([271) . classical methods for optimization problems with constraints can 
be used. We refer to [10] for a presentation of the classical Uzawa algorithm, and to [3] for 
more advanced methods. Problem ([28]) can be solved using classical methods for nonlinear 
problems, such as the Newton method (see [3]). We also refer to Chapter 7 of [18] where similar 
problems are discussed, for the classical RATTLE and SHAKE schemes used for Hamiltonian 
dynamics with constraints. 

Both schemes are consistent (the discretization error goes to when the time step At 
goes to 0) with the projected diffusion ([TSll (see [6]). Accordingly, AA„_|_i is a consistent 
discretization of //"^^ c^A^ and therefore, it can be proven [6]: 



^ T/At 

lim lim - V AA„ = F'iz) 

n=l 



which is the discrete counterpart of the trajectory average ([TQ]) . In [6], a variance reduction 
technique is proposed, which consists in extracting the bounded variation part AA^ of AA^ 
(resorting locally to reversed Brownian increments). We give some details of an adaptation of 
this method for evolving constraints in next section. 
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3.2 Discretization with evolving constraints 

When nonequihbrium dynamics are considered, the constraint is stated as (,{Qt) = z{t). The 
reaction coordinate path is first discretized as {^(0), . . . ,z{tj\fr^)} where Nt is the number of 
timesteps. For example, equal time increments can be used, in which case At = and 
tn = nAt (we refer to Remark l3.ll below for some refinements). The initial conditions Qq are 
sampled according to //So- A way to do that is to subsample a long trajectory of the projected 
SDE on So (using the schemes ([27} or ([28]) ). 

The projected SDE on evolving constraints (ETJ) is then discretized with the scheme ([21 
or ((28]) . taking into account the evolution of the constraint: 



Qn+i — Qn - VV{Qn) At + V2At/3-i Un + AA„+i Ve(g„,+i), 
where AA^+i is such that (,{Qn+i) = z(t„-|_i). 



or 



(29) 
(30) 



Qn+l — Qn - yV{Qn) At + ^/2Atp^Un + AA„+i V^(Q„), 
where AA^+i is such that ^(Qn+i) = z{tn+i)- 

It remains to extract the force part AA^_^^ from the discretized Lagrange multiplier AA„+i 
(consistently with (f22]l ). We propose two methods. First, this can be done by simply sub- 
tracting the drift and the martingale part 

Another possibility in the spirit of the variance reduction techniques used in [6] can also be 
used. Consider the following coupled dynamic with locally time-reversed constraint evolution 
(written here for the scheme (|29p ): 



Qn+l = Qn - VV{Qn) At - ^2At[j-^Un + AA^^^ V^(Q^: 



with AA^^_]^ such that: 



\ml+i)+mn+i)) = mr 



The position Q^j^i is computed as Qn+i in ([291) . but with a projection on '^2i{Qn)-^(Qn+i) 
instead of S^(f^^j), and using the Brownian increment — \/AtC/„ instead of \/AtC/„. Notice 
that in case of a constant increment for the constraints, we have i{Q^_^i) = 2^{Qn)—i{Qn+i) = 
z{tn-i)- The force part AA^^j^ is then obtained through 

AA^+i = i(AA„+i + AA^+i) (32) 

which can be shown to be a consistent time discretization of r*"+^ dkl. 

3.3 Computation of free energy using a Feynman-Kac equality 

The consistent discretization of Qt, and more precisely of //"^^ dAl, we have obtained in the 
previous section can now be used to approximate the work VV(t) defined by (f25l) by 

w„ = o, 

tn+1 '■Ji 
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using either the dynamics ((29]) or ([30]) . and the local force part of the Lagrange multiplier 
computed by ([i?T]) or ([^ . Averaging over M independent realizations (the corresponding 
works being labeled by an upper index 1 <m < M), an estimator of the free energy difference 
AF{z{T)) is, using Theorem [221 

AF{z{T)) = -r' In (j^ e-^^--) • (34) 

\ m=l / 

The estimator AF{z{T)) converges to AF{z{T)) as At ^ and M — > +oo. It is clear that 
the estimation of AF{z{T)) by (f34l) is straightforward to parallelize since the {W^^)i<m<M 
are independent. 

Notice that, even in the limit At 0, AF{z{T)) is a biased estimator. Indeed, 
exp(— /3AF(z(T))) is an unbiased estimator of exp(— /3AF(z(r))), and therefore, using the 
concavity of In, K{AF{z(T))) > AF{z{T)). Recent works propose corrections to this system- 
atic bias using asymptotic expansions in the limit M — > +oo (see for instance [24l [36] ). 

Remark 3.1 (On practical implementation). Notice that it may be useful to adaptively refine 
the time step over each stochastic trajectories, using for example the work evolution rate (Wn — 
W'n-i)n>i o,s a refinement criterion. 

As noticed in \2^^ , it is also possible to optimize the evolution of the constraint z{t), for ex- 
ample by minimizing the variance of the results obtained for a priori schedules for the evolving 
constraint on a small set of preliminary runs. 



4 Numerical results 

We present in this section some illustrations of the algorithm we have described above to com- 
pute free energy differences through nonequilibrium paths. In Section [HH a two-dimensional 
toy potential V is used, for which we can compare the results with analytical profiles. A more 
realistic test case in Section demonstrates the ability of the method to compute free energy 
profiles in presence of a free energy barrier. 

Our aim in this section is not to compare the numerical efficiency of the thermodynamic 
integration method presented in Section [1] (or any other method) with nonequilibrium com- 
putations, since it is difficult to draw general conclusions about such comparisons. However, 
we compare on a simple example in Section 14.11 the numerical efficiency of out-of-equilibrium 
computations using a few long trajectories or many short trajectories, at a fixed computational 
cost. 

4.1 A two-dimensional toy problem 

We consider the two-dimensional potential introduced in [33] 

V{x, y) = cos(27rx)(l + diy) + c^22/^ (35) 

where di and d2 are two positive constants. Some corresponding Boltzmann-Gibbs probability 
densities are depicted in Figure [H 

We want to compute the free energy difference profile between the initial state x = xq = 
—0.5 and the transition state x = xi = 0. Notice that the saddle point is (xi,yi) = (0,0) for 
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Figure 1: Plot of some probability densities corresponding to the potential (f35l) for /3 = 1, 
^2 = 27r^, and di = on the left or di = 10 on the right. 

di = 0, but is increasingly shifted toward lower values of yi as di increases. We parameterize 
the transition along the rr-axis, either with the reaction coordinate 



X — Xq 
Xi-Xo'' 



(36) 



or with the reaction coordinate (n > 2) 



Vn{x,y) 



1 + 



X — Xq 
Xi - Xo 



(37) 



For these reaction coordinates, the initial state (resp. the transition state) corresponds to a 
value of the reaction coordinate z = (resp. z = I). The analytical expression of the free 
energy difference that we consider here is, for a reaction coordinate ^{x, y) (such as ^ or -qn 
defined above) 

' r -0V{x,y)^ \ 



AF,{z) 



-p-' In 



where the distribution 6^(^x,y)-z is defined in Remark ll.ll above. Notice that even though the 
initial state Sq = {x = —0.5} and the final state T,i = {x = 0} are the same for the reaction 
coordinates ^ and r/„, the associated free energy differences differ. This is due to the fact that 
/ Vr/„, and therefore (^^(x,?/)-^ / Sr,„{x,y)-z- More precisely, 

AF^{z) = — cos(27rxo) + cos(27rxg(z)) H — — - — (cos^(27r2;o) — cos^(27rx^(z))), 

4(12 



with 



x^{z) = xo + z{xi - Xq), 
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Figure 2: Free energy profiles using the potential ([35]) with /? = 1, di = 30 and d2 = 27r^, 
and the reaction coordinate ([^^ on the left, or the reaction coordinate ([?7|l with n = 5 on 
the right. Analytical reference profiles are in dotted lines. The dashed lines (resp. the solid 
lines) represent the upper and lower bound of the 95 % confidence interval (obtained over 100 
independent realizations) for nonequilibrium computations with M = 10^ replicas (resp. with 
M = replicas). The switching time is T = 1 and the time step is At = 0.005 on the left 
and At = 0.0025 on the right. 

and 

^Fv^^i^) = -cos(27rxo) + cos(27rx^^(z)) + — — — (cos^(27rxo) - cos^(27rXr,„ (z))) 

4a2 

+ ^ln (i + ^uM^)^ 

with 

XnA^) = ^0 + ((2" - l)z + 1)1/" - l)(xi - xo). 

Free energy profiles for the two reaction coordinates considered here can then be computed 
using the discretization proposed in Section 13.31 Averaging over several realizations, error 
estimates can be proposed: in particular, the standard deviation can be computed for all 
intermediate points z € [0,1], so that, for all values z, a confidence interval around the 
empirical mean can be proposed. We represent on Figure [2] the analytical profiles, and the 
lower and upper bounds of the 95 % confidence interval for M = 10^ and M = 10^, using here 
and henceforth a linear schedule: z{t) = t/T . The initial conditions are created by subsampling 
every 100 timesteps a trajectory constrained to remain on the initial submanifold Sq. As 
announced above, the profiles obtained with t/„ and ^ are not exactly the same, though 
the general shape is preserved. These figures also show that the variance increases with z. 
Therefore, to further test the convergence of the method, it is enough here to characterize the 
convergence of the value for the end point at z = 1. 

We study the convergence of the end value AF(1) computed with the out-of-equilibrium 
dynamics with respect to the number of replicas M and the time step At, using the reaction 
coordinate ([36]) as an example. The results are presented in Table (H The time step At 
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At T M 


AF(z(T)) 


0.001 1 10^ 
0.0025 1 10^ 
0.005 1 10^ 
0.01 1 10^ 


2.056 (0.274) 
2.033 (0.259) 
2.076 (0.286) 
2.073 (0.278) 


0.005 1 10^ 
0.005 1 10"^ 
0.005 1 10^ 


2.076 (0.286) 
2.014 (0.116) 
2.001 (0.045) 



At T M 


AF{z{T)) 


0.005 1 10* 
0.005 10 10^ 
0.005 100 102 
0.005 1000 10^ 


2.014 (0.116) 
1.999 (0.029) 
2.001 (0.025) 
1.997 (0.022) 



Table 1: Free energy differences AF(1) obtained by nonequilibrium computations for the 
reaction coordinate ([^Hl) with /3 = 1, di = 1 and d2 = 30. The results are presented as follows: 

V«(™))) esU.a.. of oM.ne. 



E (^AF{z{T)) 

averages over 100 independent runs). The exact value is AF(1) = 2. 



does not seem to have any noticeable influence on the final result, as long as it remains in a 
reasonable range. As expected, the error gets smaller as M increases. 

In Tabled! we also show that, in this particular case, for a fixed computational cost and 
provided that the switching time is large enough^, computing many short trajectories is as 
efficient as computing a few longer ones (the mean and the variance are essentially unchanged) . 
This conclusion also holds for the more realistic test case presented in next section. The 
computation of many trajectories can be straightforwardly and very efficiently parallelized. 

We finally mention that we are able to exhibit the bias of the Jarzynski estimator in this 
particular case (see Section [3?3l and |36|). We observe that the estimator AF{z{T)) is generally 
greater than AF{z(T)). More precisely, averaging over 10^ realizations, with the parameters 
T = 1 and At = 0.005, we obtain the following 95 % confidence intervals for AF{z(T)), for 
various values of M: AF{z{T)) = 2.0576 ± 0.0059 for M = 10^ AF{z{T)) = 2.0095 ± 0.0026 
for M = 10^, and AF{z{T)) = 2.00075 ± 0.0010 for M = 10^ As expected, the bias goes to 
zero when M ^ oo. 



4.2 Model system for conformational changes influenced by solvation 



We consider a system composed of particles in a periodic box of side length /, interacting 
through the purely repulsive WCA pair potential [81 [28]: 



4e 
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a\6 
r 







+ e if r < ro, 
if r > ro, 



where r denotes the distance between two particles, e and a are two positive parameters and 
ro = 2^/^cr. Among these particles, two (numbered 1 and 2 in the following) are designated 
to form a dimer while the others are solvent particles. Instead of the above WCA potential, 
the interaction potential between the two particles of the dimer is a double-well potential 



Vsir) = h 



1 



(r — ro — tf)^ 



(38) 



^Of course, this threshold time depends on the system under study. 
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Figure 3: Schematic views of the system, when the dimer is in the compact state (Left), and 
in the stretched state (Right). The interaction of the particles forming the dimer is described 
by a double well potential. All the other interactions are of WCA form. 



where h and w are two positive parameters. The potential Vs exhibits two energy minima, 
one corresponding to the compact state where the length of the dimer is r = tq, and one 
corresponding to the stretched state where this length is r = ro + 2zi;. The energy barrier 
separating both states is h. Figure [3] presents a schematic view of the system. 
The reaction coordinate used is 

V-/ N \li - 92I - ro , . 

m = ^ , (39) 

where qi and q2 are the positions of the particles forming the dimer. The compact state (resp. 
the stretched state) corresponds to a value of the reaction coordinate z = (resp. z = 1). 

The parameters used for the simulations are: (3=l,e = l,a = l,h = l,w = 0.5 and = 
16. We still use a linear schedule: z{t) = t/T. The side length / of the simulation box takes 
two values: I = 1.3 (high density state) and I = 3 (low density state). Figure H] presents some 
plots of the free energy difference profiles computed using nonequilibrium dynamics, as well as 
thermodynamic integration reference profiles. The results show that nonequilibrium estimates 
are consistent with thermodynamic integration. Our experience on this particular example 
also shows that it is computationally as efficient to simulate several short nonequilibrium 
trajectories (provided the switching time is not too small, say, T ~ 1 in the units used here, 
so that the diffusion process can take place), or one single long trajectory where the switching 
is done slowly (as already observed in Section BTTI) . 

The free energy profiles highlight the relative stabilities of the two conformations of the 
dimer: at low densities (Figure HJ Left) the stretched conformation has a lower free energy 
and is thus expected to be more stable (this can indeed be verified by running long molecular 
dynamics trajectories and monitoring the time spent in each conformation). When the density 
increases, the compact conformation becomes more and more likely. At the density considered 
in Figure m (Right), the compact state already has a free energy slightly smaller than the 
stretched state. Notice also that the free energy barrier increases as the density increases, so 
that spontaneous transitions are less and less frequent. But since we know here a reaction 
coordinate, we can enforce the transition. This prevents us from running and monitoring long 
trajectories to get sufficient statistics to compare relative occurrences of both states. 
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Parameter z Parameter z 

Figure 4: Comparison of free energy difference profiles using the reaction coordinate ([Ml) , 
at low densities (/ = 3) on the left, and high densities (/ = 1.3) on the right. The double 
well potential Vs is represented in dashed line. The reference free energy difference profile 
computed with a very precise thermodynamic integration is represented in dotted line. We 
used A^TI = 101 thermodynamic integration points (uniformly distributed over (0, 1)) and 
averaged the mean force over Mti = 10^ configurations for each fixed value of z. The upper 
and lower bounds of the 95 % confidence interval (obtained over 50 independent realizations) 
for out-of-equilibrium computations are represented with solid lines. We used M = 1000 
nonequilibrium trajectories, a switching time T = 1, and a timestep At = 0.00025 (left) 
or At = 0.0005 (right). 
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A Appendix: The multi-dimensional case 

In this appendix, we generalize the previous results for nonequilibrium computation of free 
energy differences to the case of multi-dimensional reaction coordinates. 

A.l Geometric setting and basic notation and formulae. 

We consider a d-dimensional system of smooth reaction coordinates ^ = (^i, . . . ,^d) '■ 
W^, non-singular on an open domain M C 

yqeM, Tange{V^i{q),...,VCd{q))=d, 

and a smooth path of associated coordinates 

z = {zi,...,zd) : [0,T] ^M'^. 

Accordingly, we define for all t £ [0, T] a smooth submanifold of codimension d contained 
in M: 

S.(t) = G K^^, m=z{t)}cM. 

In the constraints space M.'^, coordinates are labeled by Greek letters and we use the 
summation convention on repeated indices. In the configuration space M^^, coordinates are 
labeled by Latin letters and we also use the summation convention on repeated indices. We 
denote by X -Y = XiYi the scalar product of two vector fields of M^^, by M : iV = MijNij 
the contraction of two tensor fields of M^^, and by {X (g) Y)ij = XiYj the tensor product of 
two vector fields of M.^^ . 

The d X d matrix 
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is the Gram matrix of the constraints. It is symmetric and strictly positive on ^A. We denote 
by G~l^ the (a, 7) component of , the inverse matrix of G. At each point q £ Ai, we 
define the orthogonal projection operator 

onto the normal space to "^^(q) and the orthogonal projection operator 

P = Id - 

onto the tangent space to ^^(q)- The mean curvature vector field of the submanifold is defined 
by: 

H = -V- ((detG)i/2(^-i^v^^^ {detG)-^/^VCa (40) 

and satisfies: 

We recall the divergence theorem on submanifolds: for any smooth function (j) : — > 
M?^ with compact support, 

/ divs(<^) da^^ = - I H- (/.das, (41) 

where divs(</>) = PijS/i<pj denotes the surface divergence, and as, is the induced Lebesgue 
measure on the submanifold of M^^. 

We will also use the co-area formula: for any smooth function (f> : M.^^ — > M, 

/ cl){q){detG{q))'/^dq= [ [ c^da^^dz. (42) 
These definitions and formulae are provided with more details in [6]. 

A. 2 Free energy and constrained difFusions for multi-dimensional reaction 
coordinates 

As in the one-dimensional case, the Boltzmann-Gibbs distribution restricted on the submani- 
fold T,z is defined by: 

d/i-E^ = exp{-PV)daj:,^, 

with 

Z,= [ exp{-(3V)da^^. 

JT.Z 

The associated free energy is: 

F(z) = -/?-' In (^.). 

Remark A.l (On the definition of the free energy: the multi-dimensional case). As in the 

one- dimensional case (see Remark if the particles initially evolve in a potential V , the 

classical definition of the free energy is as above, hut with V replaced by an effective potential 
V + In ((detG)^/^) . The computation of the gradient of this potential in the dynamics 
then involves second-order derivatives of ^, which can be approximated in practice by finite 
differences (see JBj). 
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For any 1 < a < d, we now introduce the local mean force along V^q, (which general- 
izes (fTTl)): 

= G~l^VC^ ■ {VV + . (43) 

As in the one-dimensional case (see Equation ([TO]) ). we obtain the derivative of the mean force 
by averaging the local mean force: 

Proposition A. 2. The derivative of the free energy F with respect to Za is given by: 

VaF{z) = f f^ d/is,. 

Proposition IA.2I is a corollary of 
Lemma A. 3. For any test function (p with compact support in M., we have: 

V„ (^j^^exp{-(3V)da^}j = ^(G-^^V^^ • - (3f^^) exp(-/3y)das.. 

Proof. It is enough to prove the formula in the case y = 0, up to a modification of the test 
function if. For any test function g : M — > M with compact support, we have (using successively 
an integration by parts on M, the co-area formula ((121), integration by parts on M"^^, and 
finally again (021)): 



/ 9{za)^a { / 'pdaj:^ \ dz = - / g'{za)(pdaj:^dz, 

= - / g °iaV (detG)^/^ dq, 



- I G-\V^^ ■ V{g o e«) ^ (detG)i/2 dq, 
I goCaV- (g-)^ if (deiGf^) dq, 

Jr3N \ J 

[ g{za) f V • (g-\V^^ if (detG)^/') (detG)-^/2 das. dz, 



which gives the result using the expression (|40]) of the mean curvature vector H. □ 
We now define the constrained diffusion (which generalizes (f2T]) ): 



(0)' 



Qo ~ /Us, 

dQt = -P{Qt)VViQt)dt + ^/2p^P{Qt)odBt + VUQt)dAf}, (44) 
dAf} = G^l^{Qt)z'^{t)dt, yi< a<d. 

The stochastic process Qt can be characterized by the following property: 

Proposition A. 4. The process Qt solution to (0ill is the only ltd process satisfying for some 
adapted ltd processes (Ai^^, . . . , A^ t)^^^ ^^] with values in W^: 

Qo ~ /^S,(o), 

dQt = -VV{Qt)dt + y^2p^dBt + VUQt)dAa,,t, 

e(Qt) = z{t). 
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Moreover, the process i-^a,t)te[o,T] can be decomposed as 
with the martingale part 



the local force part (see (J^ ) for the definition of fa) 

dAi^t = fa{Qt)dt, 

and the external forcing (or switching) term 

dA-^^ = G~^^{Qt)z'^{t)dt. 

The proof consists in computing dS,{Qt) by Ito's calculus and identifying the bounded 
variation and the martingale parts of the stochastic processes. 

A. 3 The Feynman-Kac fluctuation equality 

Theorem 12.21 is generalized as: 



Theorem A. 5 (Feynman-Kac fluctuation equality). Let us define the nonequilibrium work 
exerted on the diffusion Qt solution to (|44p by: 

f fa{Qs)z'a{s)ds= [\'^is)dAl^. 

Jo Jo 



Zz{0) 



Then, we have the following fluctuation equality: for any test function f, and\lt G [0,T], 

/ V'd/.s.(,)=IE(^(Qt)e-^^W). (45) 

In particular, we have the work fluctuation identity: \/t G [0, T], 

AF(z(t)) = F{z{t)) - F{z{0)) = -p-^ In (e (^e'^^W)) . (46) 

Proof. For any s G [0, T] and x G A^, let us introduce {Qt'^)tG[s,T]j the stochastic process 
satisfying the SDE ([441) . starting from x at time s: 



Qs — X, 

dQt'"" = -P{Qt''')VV{Qtndt+V2F^P{Qr)odBt + VUQt''')dA''al (47) 
dAf} = G-\^{Qr)z'^{t)dt, yi<a<d. 

Notice that for any s G [0,T], there is an open neighborhood {s~,s~^) x Ai^ of (5,^2(3)) in 
M X Al such that the diffusion iQt'^)te[s,T] remains in A( almost surely. This holds since this 
process satisfies d£,{Ql'^) = z'{t)dt and therefore CiQt'^) = + z{t) — z{s). This gives 
usual regularity assumptions sufficient to get a backward semi-group (t being from now on 
fixed in (0, T) and s varying in [0, t]): 

u{s, x)=E (ifiQtn exp (-(3 f fa{QrK{r) dr 
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satisfying the following partial differential equation (PDE) on (s ,s~^) x Ais- 

dsU = -Ls{u{s, .)) + f3z'^{s)faU, 
where Ls is the generator of the diffusion Qt solution to (l44l) : 

L, = : V2 - PVV ■ V + r^H ■ V + z'^{s)G-l^VCa ■ V. 

Now, using Lemma [A. 31 we have: 
d 



/ {-Ls{u{s, .)) + z'^{s)G-l,V^^ ■ Vu{s, .)) exp(-/3y)das,(,), 
- I {P~'P : VMs, •) - • Vu{s, .) + ■ Vu{s, .)) exp(-/3F)d(Ts,(,) , 

1 / (divE (Vtx(s, .) eM-py)) + ^ • V7x(s, .) exp(-/3y)) , 



= -P 
= 0, 

by the divergence theorem (|4ip . Therefore 



/ .)exp(-/3y)dcrs^,(^j = / u(0, .)exp(-/3y)(icrs^,(Q^, 



which yields 



where satisfies dSl). This proves ([^51) . and ([Ml) is obtained by taking (/? = !. □ 
A. 4 The numerical scheme 

The adaptation of the algorithm we propose for the one-dimensional case to the multi- 
dimensional case is straightforward. Indeed, the generalizations of schemes ([29]) and (f30]l 
to the multi-dimensional case are, respectively: 

Qn+l = 

where {AAa,n+i)i<a<d is such that ^(Q„+i) = z(tn+i), 



Qn+i = Qn- Vy(Q„) At + ^2Atf3~Wn + AA„,„+i Ve^(Q„), 
where (AAQ,,„+i)i<Q,<d is such that (,{Qn+i) = z{tn+i). 

The force part AA^ „ of AAq,^„ is obtained by similar procedures as those described in Sec- 
tion [321 For example, the generalization of ([3T]) is: 



The generalization of ([32]) is also straightforward. 
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Now, the estimator AF{z{T)) of the free energy difference AF{z{T)) is given by ([34ll . 
with the following approximation of the work yV{t): 

>Vo = 0, 

Wn+l - Wn^ — AA„^„_^i, 

''n+1 ''n 

which generalizes ((331) . Notice that Remark l3.ll also holds for a multi-dimensional reaction 
coordinate. 
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